#include "kppch.h"

#include "primitive.h"
#include "shape.h"

namespace kp {

	//////////////////////////////
	/// Sphere
	//////////////////////////////

	bool Sphere::Intersect(const Ray& r, float& tHit, float& epsilon, Intersection& isct, SdfRecordType& sdfRec) const {
		// Transform Ray to object space
		Ray ray = (*worldToObject)(r);
		// Compute quadratic sphere coefficients 'a, b, c'
		float a = ray.direct.dot();
		float b = 2 * ray.direct.dot(ray.origin);
		float c = ray.origin.dot() - radius * radius;
		float discriminant = b * b - 4.0f * a * c;
		if (discriminant <= .0)return false;
		// Solve high-precession quadratic equation for t values
		float rootD = std::sqrtf(discriminant);
		float q;
		if (b < 0)q = -.5f * (b - rootD);
		else q = -.5f * (b + rootD);
		float t0 = q / a;
		float t1 = c / q;
		if (t0 > t1)std::swap(t0, t1);
		if (t0 > ray.maxt || t1 < ray.mint)return false;
		float thit = t0;
		if (t0 < ray.mint) {
			thit = t1;
			if (thit > ray.maxt)return false;
		}
		// Get hit point and normal
		Point3f pHit = ray.to(thit);
		if (pHit.x == 0.f && pHit.y == 0.f)pHit.x = 1e-5f * radius;
		Normal Ng = pHit;
		if (reverseOrientation)Ng = -Ng;
		/*float cos_theta = nor.dot(r.direct);
		nor = cos_theta < 0 ? nor : -nor;*/
		// Find parametric representation of hit
		using Angle = typename float;
		Angle phi = atan2f(pHit.y, pHit.x);
		if (phi < 0.) phi += 2.f * Pi;
		Angle theta = acosf(std::clamp(pHit.z / radius, -1.f, 1.f));
		float u = phi / phi_max;
		float v = theta * 0.5f + 0.5f;
		// Set Intersection and other params
		const Transform& o2w = *objectToWorld;
		tHit = thit;
		epsilon = 5e-4f * thit;
		isct.hitp = o2w(pHit);
		isct.Ng = o2w(Ng).normalized(); // normalize!!!
		isct.Ns = isct.Ng;
		isct.uv = { u,v };
		isct.shape = this;
		return true;
	}

	bool Sphere::Sdf(const Point3f& p, float& sdfResult) const {
		return false;
	}

	BBox3f Sphere::GetBBox() const {
		return {
			{ -radius,-radius,-radius },
			{  radius, radius, radius }
		};
	}

	//////////////////////////////
	/// Triangle
	//////////////////////////////

	TriangleMesh::TriangleMesh(Ref<Transform> o2w, bool ro,
		int nt, int nv,
		const Face* vi,
		const Point3f* P,
		const Normal* N,
		const Vec3f* S,
		const Point2f* UV)
		: nTrianlges(nt), nVertices(nv), faces(vi, vi + 3 * (size_t)nTrianlges) {
		if (P) {
			p.reset(new Point3f[nVertices]);
			for (int i = 0; i < nVertices; i++)
				p[i] = (*o2w)(P[i]);
		}
		if (N) {
			n.reset(new Normal[nVertices]);
			for (int i = 0; i < nVertices; i++)
				n[i] = (*o2w)(N[i]);
		}
		if (S) {
			s.reset(new Vec3f[nVertices]);
			for (int i = 0; i < nVertices; i++)
				s[i] = (*o2w)(S[i]);
		}
		if (UV) {
			uv.reset(new Point2f[nVertices]);
			memcpy(uv.get(), UV, nVertices * sizeof(Point2f));
		}
	}

	bool Triangle::Intersect(const Ray& r, float& tHit, float& epsilon, Intersection& isct, SdfRecordType& sdfRec) const {
		Ray ray = r;
		// Calculate t values
		// o + td = (1 - b1 - b2) * pa + b1 * pb + b2 * pc = phit
		const Point3f& pa = mesh->p[mesh->faces[3 * (size_t)index].v];
		const Point3f& pb = mesh->p[mesh->faces[3 * (size_t)index + 1].v];
		const Point3f& pc = mesh->p[mesh->faces[3 * (size_t)index + 2].v];
		Vec3f e1 = pb - pa;
		Vec3f e2 = pc - pa;
		Vec3f s1 = ray.direct.cross(e2);
		float divisor = s1.dot(e1);
		if (divisor == 0.)return false;
		float invDivisor = 1.0f / divisor;
		// Calculate centroid
		Vec3f d = ray.origin - pa;
		float b1 = d.dot(s1) * invDivisor;
		if (b1 < 0.f || b1 > 1.0f)return false;
		Vec3f s2 = d.cross(e1);
		float b2 = ray.direct.dot(s2) * invDivisor;
		if (b2 < 0.f || b2 > 1.0f)return false;
		float b0 = 1.0f - b1 - b2;
		if (b0 < 0.f || b0 > 1.0f)return false;
		float thit = e2.dot(s2) * invDivisor;
		if (thit > ray.maxt || thit < ray.mint)return false;
		// Get hit point and judge if hit in triangle
		Point3f pHit = ray.to(thit);
		// Calculate normal and uv
		Normal Ng, Ns;
		Ng = Vec3f(e1).cross(e2);
		if (reverseOrientation)Ng = -Ng;
		if (nullptr != mesh->n) {
			Ns = (mesh->n[mesh->faces[3 * (size_t)index].n] * b0 +
				mesh->n[mesh->faces[3 * (size_t)index + 1].n] * b1 +
				mesh->n[mesh->faces[3 * (size_t)index + 2].n]) * b2;
		}
		else {
			Ns = Ng;
		}
		//n = Vec3f(e1).cross(e2);
		// Transmission normal?
		//float cos_theta = nor.dot(r.direct);
		//nor = cos_theta < 0 ? nor : -nor;
		Point2f uv;
		if (nullptr != mesh->uv) {
			uv = (mesh->uv[mesh->faces[3 * (size_t)index].t] * b0 +
				mesh->uv[mesh->faces[3 * (size_t)index + 1].t] * b1 +
				mesh->uv[mesh->faces[3 * (size_t)index + 2].t]) * b2;
		}
		// Set Intersection and other params
		tHit = thit;
		epsilon = 5e-4f * thit;
		isct.hitp = pHit;
		isct.Ng = Ng.normalized(); // normalize!!!
		isct.Ns = Ns.normalized(); // normalize!!!
		assert(!std::isnan(isct.Ns.x) && !std::isnan(isct.Ns.y) && !std::isnan(isct.Ns.z));
		isct.uv = uv;
		isct.shape = this;
		return true;
	}

	bool Triangle::Sdf(const Point3f& p, float& sdfResult) const {
		return false;
	}

	BBox3f Triangle::GetBBox() const {
		const Point3f& pa = mesh->p[mesh->faces[3 * (size_t)index].v];
		const Point3f& pb = mesh->p[mesh->faces[3 * (size_t)index + 1].v];
		const Point3f& pc = mesh->p[mesh->faces[3 * (size_t)index + 2].v];
		return BBox3f::Union(BBox3f(pa, pb), pc);
	}

	float Triangle::Area() const {
		const Point3f& pa = mesh->p[mesh->faces[3 * (size_t)index].v];
		const Point3f& pb = mesh->p[mesh->faces[3 * (size_t)index + 1].v];
		const Point3f& pc = mesh->p[mesh->faces[3 * (size_t)index + 2].v];
		Vec3f e1 = pb - pa;
		Vec3f e2 = pc - pa;
		return 0.5f * e1.cross(e2).length();
	}

	std::vector<Ref<Shape>> Triangle::CreateTriangleMesh(
		Ref<Transform> o2w, Ref<Transform> w2o, bool ro, 
		int nt, int nv, 
		const Face* fa, 
		const Point3f* p, 
		const Normal* n, 
		const Vec3f* s, 
		const Point2f* uv) {
		auto mesh = MakeRef<TriangleMesh>(o2w, ro,
			nt, nv, fa, p, n, s, uv);
		std::vector<Ref<Shape>> tris;
		tris.reserve(nt);
		for (int i = 0; i < nt; ++i)
			tris.push_back(MakeRef<Triangle>(o2w, w2o, ro, mesh, i));
		return tris;
	}

}